Load data

here::set_here()

## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scripts

suppressPackageStartupMessages({
  library(slingshot)
  library(tradeSeq)
  library(SingleCellExperiment)
  library(cowplot)
  library(rgl)
  library(clusterExperiment)
  library(RColorBrewer)
  library(aggregation)
  library(ggplot2)
  library(pheatmap)
  library(wesanderson)
  library(UpSetR)
  library(gridExtra)
  library(scales)
})

Import data

sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
timePoint <- readRDS("../data/finalTrajectory/timePoint_noResp_noMV.rds")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")

Chronological timepoints

plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[timePoint], alpha=.6)

df <- data.frame(UMAP1=reducedDim(sds)[,1],
                 UMAP2=reducedDim(sds)[,2],
                 time=factor(timePoint),
                 col=brewer.pal(9,'Set1')[timePoint])

ggplot(df, aes(x=UMAP1, y=UMAP2, col=time)) +
  geom_point(size=.2, alpha=.5) +
  scale_color_manual(values = brewer.pal(9,'Set1')) +
  theme_classic() +
  ggtitle("Cells colored by sampled timepoint") +
  guides(colour = guide_legend(override.aes = list(alpha = 1, size=1.5)))

# ggsave("../plots/timePointDR.png", width=9)

3D plot of trajectory

plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[timePoint], alpha=.6)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=4)

plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.6)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=4)


plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.3,
       box=FALSE)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=6)



plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.3,
       box=FALSE, axes = FALSE, xlab = "UMAP1", ylab="UMAP2", zlab="UMAP3")
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=5)
axis3d('x', at = NULL, labels = FALSE, tick = TRUE, line = 0, 
    pos = NULL, nticks = 5) 
axis3d('y', at = NULL, labels = FALSE, tick = TRUE, line = 0, 
    pos = NULL, nticks = 5) 
axis3d('z', at = NULL, labels = FALSE, tick = TRUE, line = 0, 
    pos = NULL, nticks = 5)
# rgl.postscript("../plots/tmp_trajectory.pdf","pdf")

2D plots of trajectory with cell types

# pairs(reducedDim(sds), col=brewer.pal(9,'Set1')[clDatta], pch=16, cex=1/2)

dims <- c(1,2)
plot(reducedDim(sds)[,dims], col=alpha(brewer.pal(9,'Set1')[clDatta], .6), 
     pch=16, cex=1/2, main="Dims 1 and 2")
for(ii in seq_along(slingCurves(sds))){
    c <- slingCurves(sds)[[ii]]
    lines(c$s[c$ord, dims], lwd = 3, col = "black")
}

dims <- c(1,3)
plot(reducedDim(sds)[,dims], col=alpha(brewer.pal(9,'Set1')[clDatta], .6), 
     pch=16, cex=1/2, main="Dims 1 and 3")
for(ii in seq_along(slingCurves(sds))){
    c <- slingCurves(sds)[[ii]]
    lines(c$s[c$ord, dims], lwd = 3, col = "black")
}

dims <- c(2,3)
plot(reducedDim(sds)[,dims], col=alpha(brewer.pal(9,'Set1')[clDatta], .6), 
     pch=16, cex=1/2, main="Dims 1 and 3")
for(ii in seq_along(slingCurves(sds))){
    c <- slingCurves(sds)[[ii]]
    lines(c$s[c$ord, dims], lwd = 3, col = "black")
}

Marker genes

library(ggplot2)
markers <- c("Krt5", "Krt14",
             "Sprr1a", "Trp63", #need 2 good HBC* markers here.
             "Cyp2g1", "Cyp1a2",
             "Ascl1", "Neurod1",
             "Omp")

dims <- c(1,3)
rd <- reducedDims(sds)[,dims]

markerPlotList <- list()
for(mm in 1:length(markers)){
  y <- log1p(counts[markers[mm],])
  df <- data.frame(dim1=rd[,1],
                   dim2=rd[,2],
                   y=y)
  markerPlotList[[mm]] <- ggplot(df, aes(x=dim1, y=dim2, colour=y)) +
    geom_point(alpha=.3, size=1/2) +
    scale_color_viridis_c() +
    theme_classic() +
    ggtitle(markers[mm]) +
    theme(legend.position = "none")
}


pMarker <- cowplot::plot_grid(plotlist = markerPlotList)
pMarker

Differential expression: most significantly increasing genes in each lineage

getUpGenes <- function(yhatDf, lineage){
  yhatDf1 <- yhatDf[yhatDf$lineage == lineage,]
  upGenes <- yhatDf1 %>% 
    group_by(gene) %>%
    summarize(up = yhat[20] > yhat[1]) %>%
    filter(up == TRUE)
  return(upGenes)
}

asTestRes <- associationTest(sce, lineages = TRUE)

## check which genes are increasing
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✔ tibble  3.1.0     ✔ dplyr   1.0.5
## ✔ tidyr   1.1.3     ✔ stringr 1.4.0
## ✔ readr   1.4.0     ✔ forcats 0.5.1
## ✔ purrr   0.3.4

## ── Conflicts ────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ dplyr::collapse()   masks IRanges::collapse()
## ✖ dplyr::combine()    masks gridExtra::combine(), Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count()      masks matrixStats::count()
## ✖ dplyr::desc()       masks IRanges::desc()
## ✖ purrr::discard()    masks scales::discard()
## ✖ tidyr::expand()     masks S4Vectors::expand()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::first()      masks S4Vectors::first()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce()     masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename()     masks S4Vectors::rename()
## ✖ dplyr::slice()      masks IRanges::slice()

yhatDf <- predictSmooth(sce, gene=names(sce), nPoints=20)
upGenesNeur <- getUpGenes(yhatDf, lineage=1)
upGenesSus <- getUpGenes(yhatDf, lineage=2)
upGenesHBC <- getUpGenes(yhatDf, lineage=3)

## order according to test statistic
asTestResNeurIncreas <- asTestRes[rownames(asTestRes) %in% upGenesNeur$gene,]
asTestResSusIncreas <- asTestRes[rownames(asTestRes) %in% upGenesSus$gene,]
asTestResHBCIncreas <- asTestRes[rownames(asTestRes) %in% upGenesHBC$gene,]

ooNeur <- order(asTestResNeurIncreas$waldStat_1, decreasing=TRUE)
ooSus <- order(asTestResSusIncreas$waldStat_2, decreasing=TRUE)
ooHBC <- order(asTestResHBCIncreas$waldStat_3, decreasing=TRUE)



plotlist <- list()
lineages <- c("Neur", "Sus", "HBC")
for(ll in 1:3){
  curoo <- get(paste0("oo",lineages[ll]))
  curres <- get(paste0("asTestRes",lineages[ll],"Increas"))
  plothlp <- list()
  for(kk in 1:4){
    plothlp[[kk]] <- plotSmoothers(sce, counts, gene=rownames(curres)[curoo[kk]], alpha= 1/5, lwd=3/2) +
      ggtitle(rownames(curres)[curoo[kk]]) +
      theme(legend.position = "none")
  }
  plotlist[[ll]] <- plothlp
}

pNeur <- cowplot::plot_grid(plotlist=plotlist[[1]])
pSus <- cowplot::plot_grid(plotlist=plotlist[[2]])
pHBC <- cowplot::plot_grid(plotlist=plotlist[[3]])
pNeur

pSus

pHBC

## legend
pLeg <- plotSmoothers(sce, counts, gene=rownames(curres)[curoo[kk]], alpha= 1/5, lwd=3/2) +
  scale_color_viridis_d(alpha = 1/5, labels = c("Neur", "Sus", "rHBC")) +
  guides(colour = guide_legend(override.aes = list(alpha = 1, size=1.5)))

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

pLeg

leg <- cowplot::get_legend(pLeg)

Composite main plot

p1 <- cowplot::plot_grid(NULL, pMarker, nrow=1, ncol=2, labels=c("a", "b")) #empty space for pasting 3D trajectory
p2 <- cowplot::plot_grid(pNeur, pSus, pHBC, nrow=1, labels=c("c", "d", "e"))
p3 <- cowplot::plot_grid(p2, leg, rel_widths = c(0.9,0.1))
p4 <- cowplot::plot_grid(p1, 
                         p3,
                         nrow=2)
p4

ggsave("../plots/figure1.pdf", width=12, height=10)
ggsave("../plots/figure1.png", width=12, height=10)

Shared vs unique DE genes

library(UpSetR)
asTestRes2 <- tradeSeq::associationTest(sce, l2fc = log2(1.5), lineages = TRUE)

neurGenes <- names(sce)[which(p.adjust(asTestRes2$pvalue_1, "fdr") <= 0.05)]
susGenes <- names(sce)[which(p.adjust(asTestRes2$pvalue_2, "fdr") <= 0.05)]
hbcGenes <- names(sce)[which(p.adjust(asTestRes2$pvalue_3, "fdr") <= 0.05)]
geneList <- list("Neuron" = neurGenes,
                 "Sus" = susGenes,
                 "HBC" = hbcGenes)

upset(fromList(geneList), order.by = "freq")

Figure 1 OE

Koen Van den Berge

10/1/2020

Load data

Import data

Chronological timepoints

3D plot of trajectory

2D plots of trajectory with cell types

Marker genes

Differential expression: most significantly increasing genes in each lineage

Composite main plot

Shared vs unique DE genes